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Abstract. Recently, stochastic, as opposed to deterministic, parameteriza- 
tions are being investigated to model the effects of unresolved subgrid scales 
(SGS) in large eddy simulations (LES) of geophysical flows. We analyse such 
a stochastic approach in the barotropic vorticity equation to show that (i) if 
the stochastic parameterization approximates the actual SGS stresses, then 
the solution of the stochastic LES approximates the "true" solution at appro- 
priate scale sizes; and that (ii) when the filter scale size approaches zero, the 
solution of the stochastic LES approaches the true solution. 



1. Motivation 

The immense number of degrees of freedom in large scale turbulent flows as 
encountered in the world oceans and atmosphere makes it impossible to simulate 
these flows in all their detail in the foreseeable future. On the other hand, it is 
essential to represent these flows reasonably accurately in Ocean and Atmospheric 
General Circulation Models (OGCMs and AGCMs) so as to improve the confi- 
dence in these model components of the earth system in ongoing effort to study 
climate and its variability (e.g., see [Q). Furthermore, it is very often the case that 
in highly resolved computations, a rather disproportionately large fraction of the 
computational effort is expended on the small scales (e.g., see |SJ) whereas a large 
fraction of the energy resides in the large scales (e.g., sec 6 ). ft is for these reasons 
that the ideas of Large Eddy Simulation (LES) — wherein the large scale unsteady 
motions driven by specifics of the flow are explicitly computed, but the small (and 
presumably more universal [7]) scales are modelled — are natural in this context. 

Given our interest in large scale geophysical flows with its small vertical to 
horizontal aspect ratio, we restrict ourselves to two-dimensional or quasi two- 
dimensional flows. Previous models of the small scales in how they affect the 
large scales in the momentum equations or equivalently the vorticity equation in 
incompressible settings have mostly been confined to an enhanced eddy viscosity 
or nonlinear eddy viscosity like that of Smagorinsky or biharmonic viscosity (e.g., 
see [HI El EH EH EH El)- Given the non- unique nature of the small scales with re- 
spect to the large scales [5], the aforementioned use of deterministic and dissipative 
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closures seem rather highly restrictive. On the other hand, it would seem desirable 
to actually represent a population of eddies that satisfy overall constraints of the 
flow rather than make flow specific parametric assumptions. This has led to recent 
investigations of the possibility of using stochastic processes to model the effects 
of unresolved scales in geophysical flows (e.g., see More recently, subgrid 

scale (SGS) stresses have been analysed in simple but resolved flows as a possible 
way to suggest stochastic parameterizations as e.g, in 0] 1301 EH] . These efforts 
have been preceded, of course, by various attempts to model anomalies in geophys- 
ical flow systems as linear Langevin equations (e.g, |281 121 ITrjl \26\ ) and the analysis 
of stochastic models in isotropic and homogeneous three dimensional turbulence 

(e.g., pang). 

In this paper, we analyze the stochastic approach to parameterization in the 
barotropic vorticity equation and show that (i) if the stochastic parameterization 
approximates the SGS stresses, then the stochastic large eddy solution approxi- 
mates the "true" solution at appropriate scale sizes; and that (ii) when the filter 
scale size approaches zero, then the solution of the stochastic LES approaches the 
true solution. In the next section we present a set of computations that demon- 
strates the use of stochastic parameterizations and in §3, we prove the main results 
on approximation and convergence of LES solutions using stochastic parameteriza- 
tions. 

2. Stochastic parameterization 

We consider the simple setting of the beta-plane barotropic vorticity equation 
(equivalently the two-dimensional (2D) quasi-geostrophic (QG) model) |271 I34| : 

(2.1) q t + J(4>, q) + ffl)x = f(x, V, t) + vAq - rq, 

on a bounded domain D with piecewise smooth boundary 3D. Here the vorticity 
q(x,y,t) is given in terms of streamfunction ip(x,y,t) by q — A-0. (3 is the merid- 
ional gradient of the Coriolis parameter, v > the viscous dissipation constant, r 
> the Ekman dissipation constant and f(x,y,t) the wind forcing. The forcing / 
is always assumed to be mean-square integrable both in time and in space. In addi- 
tion, A = d xx + d yy is the Laplacian operator in the plane and J{h, g) = h x g y — h v g x 
is the Jacobian operator. The boundary condition (BC) is q — 0, ip — on dD 
and initial condition (IC) is q(x 1 y,0) = qa(x,y). 

Fine mesh simulations (q) are used to obtain the benchmark solution q through 
convolution with a spatial filter G$(x, y), with spatial scale 6 > 0: 

(2.2) q(x,y,t) :=q*G S 

(a; 2 +y 2 ) 

We use a Gaussian filter Gs(x,y) = -^e ^ , where S > is the filter 
size and the filter is such that (1) q * G$ is infinitely differentiable in space and 
(2) q * Gs — * q as S —> in L 2 (D). Note that the Fourier transform of Gs is 

^ 5 2 (fc 2 +fc§) ^ ^ 

Gs(ki,k 2 ) = e 4 , and that q * G s (ki,k 2 , t) = Gs(ki,k2)q(ki, k 2 ,t). 
On convolving (|2.1|l with Gs the large eddy solution q is seen to satisfy 

(2.3) q t + J$, q) + 0$ x = f(x, y, t) + isAq - rq + R(q), 
where the SGS stress term R(q) is defined as 



(2.4) 



R(q) :=J$,q)-J(il),q). 
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Note that since q ^ q, the SGS stress term R(q) above is more than the Reynolds 
stress and the SGS stress is usually further divided into three components, the 
explicit Leonard stress, and the cross stress and the SGS Reynolds stress that 
require further modeling. However, for our purposes, we will consider R(q) in its 
entirety. Since R(q) depends on q as well as q, the equation l|2.3|> is not a closed 
system. We need to model or prescribe R{q) in terms of resolved quantity q. On 
the other hand, R(q) may be explicitly diagnosed from a fully-resolved run, given 
the filter. 

An analysis of R(q) reveals that its time-mean is much smaller than its stan- 
dard deviation, and that its temporal behavior is highly irregular, leading to the 
possibility of approximating R(q) by a suitable stochastic process a(q, u>) (defined 
on a probability space (fl,J~, P), with u £ Cl, the sample space, a— field T and 
probability measure P). With such a putative stochastic closure, the LES model 
becomes a random partial differential equation (PDE) for Q ~ q: 

(2.5) Q t + J(% Q) + && x = />, y, t) + vAQ - rQ + a(Q, to), 

where Q = A* and f(x,y,t) := / * G s , with BC: Q = 0, * = on dD and 
IC: Q(x,y,0) = q~o(x,y). Note that when the stochastic process a(Q,u>) does not 
depend explicitly on Q — as is the case when for example, either R(q) or some of 
its statistical properties are used to construct <j(ui) — we end up with an additive 
stochastic closure. The more general case of the stochastic closure wherein there is 
an explicit dependence on the state of the system Q, corresponds to a multiplicative 
stochastic closure. 

While one can get an idea of the stochastic forcing to be used to represent the 
effects of unresolved subgrid scales in the LES runs by analysing resolved runs at 
the scale of the LES computation, the selection of a specific functional form for the 
stochastic parameterization is beyond the scope of the present article. We aim at 
providing a quantitative guide to selecting the stochastic parameterization. 

2.1. Numerical experiments. We now briefly present a set of computations us- 
ing the beta-plane barotropic vorticity equation (|2.1|l in a rectangular midlatitude 
basin. Finite differencing in space is used along with Runge-Kutta time stepping, 
and other details of the setup may be found elsewhere 22- The steady forcing is 
uniform in the zonal (x) directions and sinusoidal in the meridional (y) direction 
corresponding to a double-gyre wind forcing. At the parameter values that we are 
presently consider, the circulation is highly variable, but statistically stationary. 
We therefore consider long time averages over the attractor in place of ensemble 
averaging (over u>). 

Fig. 1 shows the contour plots of the time average of streamfunction and po- 
tential vorticity as they emerge in the resolved computations. A discussion of the 
phenomenology of this circulation may be found in |82| . This simulation is then 
analysed using a Gaussian filter with a width that is four times the grid spacing of 
the resolved computations to obtain the SGS stress R(q). 

Next we consider a pair of coarse-scale simulations in which the grid spacing 
is four times that of the resolved computation in both directions. Fig. 2 shows 
the time average of the streamfunction and potential vorticity as emerges from the 
coarse-scale computation when a(Q, to) in H2.5(l is set to zero. The main differences, 
as compared to the resolved runs, clearly are the absence of the outer gyres in the 




Figure 1. Resolved simulation. Contour plots of the average 
over the attractor (time-average) of streamfunction and potential- 
vorticity. Contour intervals are the same in all the figures and are 
indicated on top of the figures in the form (min, max, increment). 



streamfunction field and the large amplitude grid scale oscillation in the potential 
vorticity field. 

In the next case, we set the statistics of tr((5,w) (viz., its spatial and temporal 
correlation functions, amplitude and probability distribution function) identical to 
those of R(q) previously diagnosed from the resolved simulation. This we do by 
using the actual time history of the SGS stress R(q) in a coarse-resolution run 
in which the initial condition is slightly perturbed from the initial condition of 
the resolved run from which R(q) is diagnosed. Given the highly chaotic nature 
of the flow, the effect of the initial perturbation is to quickly lead to a complete 
decorrelation of the SGS stress forcing term R(q) from the state of the system Q. 
Thus, the actual time history of the diagnosed SGS stress R(q) supplied to the 
coarse-resolution run acts effectively as an additive stochastic closure a(uj) of the 
SGS stresses in this LES. The time averages of the circulation for this case is shown 
in Fig. 3. Comparing Fig. 2 and Fig. 1 and Fig. 3 and Fig. 1, it is clear that the 
differences between the resolved case and the case with a(Q 7 w) similar to R(q) 
is much smaller than the differences between the resolved case and the case with 
a(Q,uj) = 0. 

These numerical experiments suggest that stochastic parameterizations can pro- 
vide good representations of subgrid scales. So, the question then is as to how such 
an approach can be justified. 
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Figure 2 . Coarse unresolved simulation. Contour plots of stream- 
function and potential- vorticity, as in Fig. 1. The stochastic pa- 
rameterization term a(Q,io) is set to 0. 



In an attempt to answer this question, we will prove, following the approaches 
in |191 I29j in our stochastic context, that, under appropriate conditions on the 
stochastic parameterization (that appear easier to check in our case, c.f., Theorem 
1 part (i) below and the Assumption A2 in [29 ): 

(2.6) E\\q- Q\\ 2 < C(u,r,q ,T) ■ E / \\R(q) - a(Q, ^)|| 2 , < t < T, 

Jo 

(2.7) E||g-Q]| 2 0, as S -> 0, < t < T. 

where C(-) > is a constant, [0, T] is the computational time interval, ¥.(Z(ui)) := 
J n Z(oj)dF(u)), and || • || is the (spatial) norm in the space L 2 (D) of spatially mean- 
square integrable functions: L 2 (D) := {/ : ||/|| = J J D f(x,y)dxdy < oo}. 



3. Main results 

Standard abbreviations L 2 = L 2 (D), Hq = Hq(D), k = 1, 2, . . ., are used for the 
common Sobolev spaces in fluid mechanics |3fil I22j . with < •, • > and || • || denoting 
the usual (spatial) scalar product and norm, respectively, in L 2 {D): 



<f,9>-=J fgdxdy, ll/H := \/< f,f > = J J f(x,y)dxdy 
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Figure 3. LES with stochastic parameterization. Contour plots 
of streamfunction and potential- vorticity, as in Fig. 1. The statis- 
tics of a(Q,ui) are identical to those of R(q). 



We need the following properties and estimates (see also [5311201) of the Jacobian 
operator J : H$ x Hq — > L 1 : 

/ J{f,g)hdxdy = - J(f,h)gdxdy, J(f,g)gdxdy = 0, and 
Jd Jd Jd 



J(f,g)dxdy 



D 



< ||V/||[|Vg|| for all f,g,h£ i?g, 



J(Af,g)Ahdxdy 



< 



2\D\ 



||A/|| ||A 5 || \\Ah\\ for all f,g,he H 2 . 



In addition, we recall the Poincare, Young, and Gronwall inequalities below: 

||<7]] 2 - f g 2 (x,y)dxdy < ^ f \Vg\ 2 dxdy = ^|]V.g|| 2 for g e Hi 
Jd tt J d tt 

AB < -A 2 + —B 2 for A, B > and e > 0. 
~ 2 2e 

Assuming that y(t) > 0, <7(i) and h(t) are integrable, if || < gy + h ioi t > to, then 



y(t) <y(to)e ! '° 9{T)dT + [ h(s)e~^ 9{T)dT ds, t>t . 
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Lemma 3.1. ( 20 ) The quasi- geostrophic motion described by 12.1)) satisfies the 
enstrophy estimate : 

(3.1) |M| 2 < \\ qo \\ 2 e- 2at + - f \\f(s)\\ 2 e 2a ^ds, < t < co 

r Jo 

(3.2) < \\q \\ 2 e 2 ^ T + -e 2 ^ T T \\f(s)\\ 2 ds, 0<t<T 

r Jo 

where a = | + p=jj — ||/3| (^^- + 1^ luii/i |D| denoting the area of the domain D. 
Note that a is positive in the case of no rotation ([3 — 0). 

Theorem 3.2. (i) Stochastic Approximation: 

If the stochastic paramerization a(Q,u>) is such that 

(3.3) / \\o~(Q,u>)\\ 2 dt < A1(T), almost surely for oj G ft, 
Jo 

for some constant M > depending on computational time interval, then 

(3.4) E||g-Q|| 2 < C(v,r,qo,T) ■ E ( \\R(q) - a(Q,uj)\\ 2 dt, < t < T, 

Jo 

for any fixed time interval < t < T. 

This implies that, if the stochastic paramerization a(Q,ui) approximates the SGS 
stress R{q), then the LES solution Q approximates q, in mean-square sense, 
(ii) Scale convergence: If the stochastic paramerization a(Q,uj) satisfies 

(3.5) E / \\a{Q,uo)\\ 2 dt -> 0, as 5 -> 0, 

Jo 

for all LES solutions Q of \2.5\) . then 

(3.6) E||g-Q|| 2 -> 0, as 5 -> 0, < t < T. 

This implies that, if the stochastic paramerization u{Q,uj) becomes smaller (collec- 
tively in computational time interval) as the cut-off scale size 8 decreases, then the 
LES solution Q approximates the original solution q better, in mean-square sense. 

Remark 3.3. Condition l|3.3|) means that the stochastic paramerization a{Q,ui) is 
square-integrable in time and space, and its norm in the space L 2 ((0, T); L 2 (D)) is 
almost surely bounded on the computational interval. 

Remark 3.4. Condition (|3.5|l means that the variance of the stochastic parameriza- 
tion a(q,u>), collectively in the finite time interval of numerical simulation, becomes 
smaller and smaller as the cut-off scale size <5 decreases. 

To prove part (i), denote U = q-Q, so that U = A(V>-*). Note that f7(0) = 0. 
Subtracting (|2.3|) from (|2.5p . we see that U satisfies 

(3.7) U t = -J{i>, q) + J{% Q) - 0$ x - + isAU - rU + [R(q) - a(Q, w)]. 

Multiplying this equation by U , integrating over D and noting that q = U + Q, we 
obtain 

(3.8) = -J J ^ ~ *' ® U ~ Pf ~ **) U 

(3.9) -H|VC/|| 2 -r||C/|| 2 + f [R{q)-a(Q,u)]U. 

JD 
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Note that, for < t < T, 

J(ip -y 7 Q)Udxdy 



D 



< 



2\D\ 



WQW Wf 



(3.10) < 



2\D\ 



{\\qo\\ 2 e 



2 2\a\T _ 



(\\m\\ 2 + MQ^)\\ 2 )d s } \\u\\< 



where we used the Lemma 13.11 on the LES model l|2.5|) . Also, by the Young and 
Poincare inequalities |36| we have 



(i (ip x - ^ x )Udxdy 



D 



< 



< 



(ip x - dxdy + U dxdy 



D 
7T 



D 



U 2 dxdy+ / U 2 dxdy 



/? / (<tp x -V x )Udxdy 



< 



1 



that is 
(3.11) 
Moreover 
(3.12) 

Putting all these estimates into 1)3. 8[l . we obtain 
d „„„, / I I2\D\ 



n 



D 



i WW 



[R(q) - a(Q,u)]U\ < h\R(q) a(Q,uj)\\ 2 + \\\U\\ 2 . 

D 1 z 



5^ 2 \- a+1 2 



(3.13) 



\M 2 e 



2 2\ a \T 



1 



,2|q|T 



+ ||fl(< ? )-a(Q,a/)|| 2 



Notice that a is defined in Lemma 13.11 in terms of physical parameters. By the 
Gronwall inequality |3fi| and noting U(0) = 0, we obtain 

(3.14) E||g-Q|| 2 = E\\U\\ 2 < C(v,r,q ,T)- E f \\R(q)^a(Q,uj)\\ 2 dt, < t < T, 

Jo 

where C > is a constant. This proves part (i) of Theorem 1. 

To prove part (ii) of Theorem 1, denote V = q — Q, so that V — A(r/> — "J). 
Subtracting equation l|2.1[l from l|2.5(l leads to 

(3.15) V t = -J(ip, q) + J(% Q) - P(ip x - * B ) + i/Ay - rV + (/ - /) - a(Q, u). 
Similar to the approach in proving part (i) above, we estimate 

fT \ 



2\D\ 



ho\\ 2 e 



2 P 2 M T _i_ i p 2|a|T 



ll/WII^S 



(3.16) 



+ ||/-/|| 2 + HQ,c)|| 2 
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By the Gronwall inequality again, we obtain 

n\Vf < C 1 ( l /,r, go ,T)E|| (Z o-go|| 2 
(3.17) + C2(v,r,q ,T)- E f [\\f - f\\ 2 + \\a(Q,u>)\\ 2 ] dt, < t < T, 

JQ 

where C\, C2 are positive constants. Due to the property of Gs, both \\qo — q~o\\ and 
||/ — /|| go to zero as S — > 0. Together with the condition lj3.5JI . we finally see that 
mV\\ 2 = E\\q - Q\\ 2 -> as 5 -> 0, completing the proof of Theorem 1. 

We view Theorem 1 as a starting point in our study of stochastic parameter- 
ization of SGS stress related terms in the LES of geophysical flows in the sense 
that it comments on the dependence of the LES solution only on the variance 
of the stochastic parameterization (its closeness to that of the actual SGS stress 
terms). Thus, it may be argued that the stochastic nature of the parameterization 
is not central to the results of this paper. However, to the extent that variance is 
one of the most fundamental characteristics of a stochastic process, understanding 
the dependence of the LES solution on it is important. What is now desirable is 
further characterization of the dependence of the LES solution on various other as- 
pects of the stochastic parameterization such as its temporal and spatial correlation 
structure and its probability distribution function. These are subjects of ongoing 
research and we hope to report on them in the future. 
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